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MULTIGRID METHODS FOR HELLAN-HERRMANN-JOHNSON 
MIXED METHOD OF KIRCHHOFF PLATE BENDING PROBLEMS 

LONG CHEN*, JUN HUt, AND XUEHAI HUANG* 


Abstract. A V-cycle multigrid method for the Hellan-Herrmann-Johnson (HHJ) discretization 
of the Kirchhoff plate bending problems is developed in this paper. It is shown that the contraction 
number of the V-cycle multigrid HHJ method is bounded away from one uniformly with respect to 
the mesh-size. The key is a stable decomposition of the kernel space which is derived from an exact 
sequence of the HHJ method. The uniform convergence is achieved for V-cycle multigrid method 
with only one smoothing step and without full elliptic regularity. Some numerical experiments are 
provided to confirm the proposed V-cycle multigrid method. The exact sequences of the HHJ method 
and the corresponding commutative diagram is of some interest independent of the current context. 


Key words. Kirchhoff plates, Hellan-Herrmann-Johnson mixed method, multigrid method, 
exact sequence, stable decomposition 


1. Introduction. We consider multigrid methods for solving the saddle point 
system arising from the Hellan-Herrmann-Johnson (HHJ) mixed method discretiza¬ 
tion (cf. [28, 29, 31]) of a fourth order equation: the Kirchhoff plate bending problem. 

Linear systems arising from discretization of fourth order partial differential equa¬ 
tions are difficult to solve due to the poor spectral properties. For conforming 
finite element methods of the biharmonic equation, some multigrid methods are stud¬ 
ied in [50, 10, 48, 52]. In practice since it is hard to construct finite elements, 
nonconforming finite element methods (cf. [20, 33, 44]), notably the Morley element 
(cf. [35, 45, 44, 46]), Zienkiewicz element (cf. [6, 43]) and Adini element (cf. [1, 33, 44]), 
are favored for the fourth order equation. Optimal-order nonconforming multigrid 
methods with the full regularity assumption are developed in [11, 37, 53, 49, 39]. 
Without assuming full elliptic regularity, similar results are obtained in [42, 13, 51]. 
For (7° interior penalty methods of fourth order equations in [14, 23], it is proved 
in [15] that V-cycle, F-cycle and W-cycle multigrid algorithms are uniform contrac¬ 
tions. Standard mutligrid solvers for the Poisson operator are used to design effi¬ 
cient smoothers. An algebraic multigrid method by smooth aggregation is developed 
for the fourth order elliptic problems in [41]. In all of these work, special intergrid 
transfer operators are necessary for both these conforming and nonconforming multi¬ 
grid methods, since either the underlying finite element spaces are non-nested or the 
quadratic forms are non-inherited. The contraction number of V-cycle, W-cycle or 
F-cycle multigrid method can be proved to be less than one uniformly with respect 
to the mesh level provided that the number of smoothing steps is large enough. 

We shall develop a multigrid method for the Hellan-Herrmann-Johnson discretiza¬ 
tion of the Kirchhoff plate bending problems in the mixed form. The resulting linear 
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system is in the following saddle point form 


( 1 . 1 ) 


(M 




which is considered harder to solve than the symmetric positive counterpart due to the 
indefiniteness of the saddle point system. To this end, the hybridization technique 
is applied to the HHJ mixed method by introducing a Lagrange multiplier, which 
changes the saddle point system to a symmetric positive definite (SPD) problem 
(cf. [26, 2]). It is shown in [2] that the resulted SPD problem in the lowest order HHJ 
method is equivalent to a modified Morley method. As we mentioned earlier, however, 
multigrid algorithms for Morley element method have been only proved to be optimal 
with special intergrid transfer operators and large enough number of smoothing steps. 

We shall apply the approach developed in [19] to design an effective multigrid 
methods for solving (1.1). The smoother of our multigrid method is a multiplicative 
Schwarz smoother based on a multilevel decomposition of the null space ker(iJ). Since 
the finite element spaces of the HHJ method are nested, the coarse-to-fine intergrid 
transfer operator are simply the natural injection. 

The key to the analysis and the algorithm is a stable multilevel decomposition 
of the null space ker(divdiv). To this end, we first establish exact sequences for the 
HHJ mixed method of Kirchhoff plates in both continuous and discrete levels. Af¬ 
ter achieving a decomposition of the finite element space for the stress based on the 
discrete exact sequence for the HHJ method, a stable decomposition and the strength¬ 
ened Cauchy Schwarz inequality are derived using the standard technique as in [47]. 
Then according to the theoretical results developed in [19], the contraction number 
of our V-cycle multigrid HHJ method is bounded away from one uniformly with re¬ 
spect to meshsize with even only one smoothing step. Since a stable decomposition 
is obtained using L^-projection, the full regularity assumption is not needed neither 
in our approach. As far as we know, our V-cycle multigrid method is the first work 
possessing these two merits among the multigrid methods for solving the fourth order 
partial differential equation directly. 

Although the multigrid method used here and its convergence follow from the 
framework developed in [19], this example has special feature which lead to rather 
difficult analysis than examples listed in [19]. Furthermore, the Hilbert complex for 
HHJ method revealed in this paper is of some interest independent of the current 
context and will play a central role in the design and analysis of the HHJ method [3] , 
c.f. the convergence of adaptive finite element methods for HHJ method established 
in [30]. We emphasize this contribution by listing the commutative diagram for the 
HHJ method as follows. Details on the spaces and interpolation operators can be 
found in Section 2.2. 


Pi(D;M2)^ 


H\n- 


V^x 


ff^(divdiv,D;S) 


divdiv 


iJ-i(D) 


Pi(D;] 


■5, 


Ih 


V^x , 

(divdiv);. J 


Vh 


■n 


The rest of this paper is organized as follows. The HHJ method for Kirchhoff 
plates and the corresponding exact sequence and commutative diagram are present 
in Section 2. Then we construct stable decomposition and prove the strengthened 
Cauchy Schwarz inequality for HHJ method in Section 3. In Section 4, we show and 
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analyze the V-cycle multigrid method for HHJ method. Some numerical experiments 
are given to testify our multigrid method in Section 4 as well. 

2. Mixed Methods for the Plate Bending Problem. Assume a thin plate 
occupies a bounded simply connected polygonal domain C Then the mathe¬ 
matical model describing the deflection u of the plate is governed by (cf. [25, 38]) 

{ Ccr = —in O, 
div diver = — / in id, 
u = dnU = 0 on did, 

where n is the unit outward normal to did, V is the usual gradient operator, div div 
stands for the divergence operator acting on vector-valued (tensor-valued) functions 
(cf. [38]). Here, C is a symmetric and positive definite operator defined as follows: for 
any second-order tensor r, 

1 v , 

Ct ■= - -T - --^(tr-r Z 

1 — V 1 — 

with I a second order identity tensor, tr the trace operator acting on second order 
tensors, and v G L°°(id) the Poisson ratio satisfying inf !^ > 0 and supz/ < 0.5. 

2.1. Hellan-Herrmann-Johnson Method. Denote the space of all symmetric 
2x2 tensor by S. Given a bounded domain G C and a non-negative integer m, 
let H'^{G) be the usual Sobolev space of functions on G, and X) be the usual 

Sobolev space of functions taking values in the finite-dimensional vector space X for 
X being S or . The corresponding norm and semi-norm are denoted respectively by 
j| • ||m,G and j • \m.G- If G is D, we abbreviate them by jj • \\m and j • jm, respectively. Let 
H^{G) be the closure of C(j“(G) with respect to the norm j| • llm,G- Pm{G) stands for 
the set of all polynomials in G with the total degree no more than m, and Pm{G;'K) 
denotes the tensor or vector version of Pm{G) for X being S or respectively. 

Let {Th}h>o be a regular family of triangulations of D. For each K G Th, denote 
by tik = ini,n 2 )'^ the unit outward normal to dK and write tic := (^ 1 ,^ 2 )^ = 
(—n 2 ,ni)^, a unit vector tangent to dK. Without causing any confusion, we will 
abbreviate and tx as n and t respectively for simplicity. Let £h be the union 
of all edges of the triangulation Th and the union of all interior edges of the 
triangulation Th- Set for each K G Th 

£h{K) -{eGSh-eC dK}, £1{K) := {e G e C dK}. 

For any e G Eh-, hx a unit normal vector rig := {ni,n 2 )’^ and a unit tangent vector 
te '.= (—n 2 ,ni)^. For a column vector function cj) = (^ 1 ,^ 2 )^) differential operators 
for scalar functions will be applied row-wise to produce a matrix function. Similarly 
for a matrix function, differential operators for vector functions are applied row-wise. 
Discrete differential operator div^ is defined as the elementwise counterpart of div 
with respect to the triangulation Th- For a second order tensor-valued function r, set 

Mn{T) := rTTn, M„t(x) := t^rn, 

on each edge e G Eh- Next, we introduce jumps on edges. Consider two adjacent 
triangles K^ and K~ sharing an interior edge e. Denote by and n~ the unit 
outward normals to the common edge e of the triangles K'^ and K~, respectively. 
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For a scalar-valued function w, write u"*" := v\k+ and v := v\k- ■ Then define jumps 
on e as follows: 


[u] := v^Ue ■ n'^ + v~ne ■ n~. 

On an edge e G K lying on the boundary 90, the above term is defined by 

[u] := vrie ■ riK- 


For any second order tensor-valued functions cr and t, set 

2 

(T : T := CTijTij. 




Then we define some Hilbert spaces. Based on the triangulation Th, let 

S:={4>€ : c/)|k e ViF € %} , 

V := {r e : t\k e H\K-,S) yK e % and [M„(T)]|e = 0 Ve S 

V:={v€H^{n):v\KyH^iK) WK £ %} ■ 

The corresponding finite element spaces are given by 

Sh := {0 G : 01^ G Pr{K-,m.^) VK G %} , 

Vh:={T£V-.T\K£Pr-i{K-,S) yxeTh}, 

Vh ■■= {v £ H^in) : v\k £ PriK) WK £ %} 
with integer r > 1. 

With previous preparation, the Hellan-Herrmann-Johnson (HHJ) mixed method 
(cf. [28, 29, 31]) for problem (2.1) is given as follows: Find {cTh,Uh) £ Vh Ph such 
that 

( 2 . 2 ) aicrh,T)+ biT,Uh) = 0 VT£Vh, 

(2.3) b{(Th,v) =-[ fvdx yv£Vh, 

Jn 

where 

a{cr,T) := / Ccr:Tdx V(t,tGV, 

Jn 

b{T,v):=— / (div/ix) • Vudcc -I- ^ / Mnt{’T)dtvds VxGV,uG7^. 

Jo. KeTh 

The boundary condition for the deflection u = 0 on dfl is imposed into the space 
Ph whereas the boundary condition for the rotation dnU = 0 on 9H is imposed weakly 
in the variational form (2.2). If the plate is simply supported along the boundary, i.e. 
the boundary condition is now u = 0, M„(cr) = 0 on dfl, we only need to modify V as 

Vo := {x G V : Mn{T) = 0 on 9H} . 


It was shown in [5, 24, 8] that the HHJ method (2.2)-(2.3) is well posed. And the 
inf-sup condition holds as follows (cf. [30, Lemma 4.2]) 


\\Vh\\2,h < sup 


b{'^h,Vh) 

\\Th\\o,h 
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y Vh£ Ph- 




where mesh dependent norms are 


K^Th 

hWl.h ■■= Ikllo + X ^e||M„(T)||g 

eGffe 

And it possesses the optimal a priori error estimates provided cr and u are smooth 
enough: 


Ik - CT/illo ^ ^’'Ikllr, 

Ik - ^i/illl k ^’'(Ikllr- + Ikllr+l)- 

Reliable and efficient a posteriori error estimators, as well as the convergence of 
an adaptive HHJ method, can be found in [30]. 

2.2. Hilbert Complex for HHJ Methods. In this section, we shall derive 
the exact sequence and commutative diagram for the HHJ method (2.2)-(2.3). 

For a vector-valued function cf) = {(l)i, 4>2)'^, denote by (f>'^ := {—4>2,4>i)'^ the 
vector perpendicular to (p. The standard symmetric gradient operator is 

£(<A) = ^(V<A+(V0f). 

The symmetric curl operator will be defined analogically 

V® X </) := - (curl(/) -1- (curl0)^) . 

Let 

It is easy to see that is exactly the rigid body motion space where 

:= {k e £ Pi(H;Rk}- 

Lemma 2.1. The following sequence for Kirchhoff plates 


_ C V^x divdiv 

(2.4) Pi(H;Rk-->0 


is an exact complex. 

Proof. By direct computation, it is easy to see that (2.4) is a complex, i.e. 
V® X (Pi) = 0 and divdiv V®x =0. We then verify the exactness. 

Let us first show that ker(V®x) = Pi(H;M^). For any 4> G C°°(r2;K^) satisfying 
V® X 0 = 0, it holds 

V® X (^ = L^e(4,^)L = 0. 
where X = ^ ^ 0^ ) Thus we have 

e{cp^) = 0, 
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which implies (f> G Pi(r2;K^). 

Next we demonstrate that ker(divdiv) = V® x using the similar argu¬ 
ment adopted in [7, Lemma 1] and [30, Lemma 3.1]. First of all, V® x C 

ker(divdiv) by direct computation. For any r € ker(divdiv), there exists v G 
such that divx = curlu = —div(uL), which implies div(T -|- vL) = 0. Hence there 
exists a vector function 4> G satisfying 

T + vL — curl0. 

Since r is symmetric, we have r = V® x (/>. Thus ker(divdiv) C V® x 

Finally we show that divdivC“(H; S) = By the elasticity complex in [4, 

p. 405], the divergence operator div : C°°(H;S) —is surjective. And due 
to the de Rham complex in [3, p. 27], the divergence operator div : —)• 

C°°(r2) is also surjective. Hence we have divdivC°°(H; S) = C°°(H). □ 

We then derive an exact sequence with less smoothness. To this end, we define 
as 


{Bt,v) := b{T,v) 'ivGV. 

For any (x,u) G V x V with v G it follows from integration by parts and the 

fact [M„(t)] jgi = 0 that 

{Bt,v) = ( T :V^vdx — ^ f Mn{T)dnvds 

Jn KeTh 

(2.5) ^ (divdivT,u)^-2(o)xff2(f2). 

On the other side, for any (t,u) GV xV with t G iT(div, H;S) := {r G L^(H;S) : 

divT G it follows from the fact [M„t(T)] = 0 that 

{Bt,v)=— / (divx) • Vudcc-b ^ / Mnti'r)dtv ds 
Jn KGTh 

= - = (divdivT,w)^-i(f2)xj^i(f2). 

Therefore the bilinear form b{-,-) can be defined either on ff(div,r2;S) x Hq{^) as 
B = div div in distribution sense or L^{Vl) x Hq{Q) with B — div div in 

distribution sense. However, conforming finite element spaces of H (div, H; S) 
or i7Q(H) are difficult to construct. We strike a balance of the smoothness of these 
two spaces and understand the bilinear form 6(-, •) being defined on V x P and thus 

divdiv : ff“^(divdiv, H; S) —> 

with space iT“^(divdiv, H;S) := {t G L^(H;S) : divdivr G which was 

firstly introduced in [36]. 

Making use of the similar argument as in Lemma 2.1, we can acquire a Hilbert 
sequence for Kirchhoff plates as follows. 

Lemma 2.2. The following Hilbert sequence for Kirchhoff plates 

( 2 . 6 ) 

_ C V®x divdiv 

Pi (H; K^)- {Gi- m 2)-►Ff-i(divdiv, H; S)- 
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is an exact complex. 

Remark 2.3. A less smooth exact Hilbert sequence for Kirchhoff plates is 
(2.7) 

_ C V'^x divdiv 

- ^H-^{dwdW,n-S) -->0, 

where i?“^(divdiv,H;S) := {r G : divdivr G Finite element 

spaces of Ff“^(divdiv,H;S) is, however, difficult to construct. Indeed in the HHJ 
method, the space V and Vh are not subspaces of ff”^(divdiv, H; S) neither. That 
is, HHJ method is still a non-conforming method. □ 

Remark 2.4. The dual complex of (2.7) is 

C rot 

0-- ^Hoirot, H; S)- >Ll{n; -^q. 


where 

JTo(rot, H; S) := {t G T^(H; S) : rotr G T^(H; R^), and rt = 0 on 9H}, 


Tg(H;R2) := {(/>G i^(H;R2) : [ (f)dx = 0}. 

Jq 

It is interesting to notice that the last exact sequence is an rotation of the elasticity 
complex in two dimensions [4, (2.1)]. □ 

In the discrete level, we shall derive a similar exact sequence for the finite element 
spaces introduced before. To this end, we first discuss the discretization of the two 
differential operators x and divdiv. Since V® x only requires the smoothness, 
it can be naturally discretized by choosing the finite element space C iJb The 
difficulty is the discretization of operator divdiv. First we can understand B : Vh ^ 

K as 

{Bt,v) := b{T,v) y V G Vh- 

Using the Riesz representation induced by the L^-inner product, we can identify 
with Vh and finally define (divdiv)/i ■ Vh Vh as follows: for any r G Vh^ 
(divdiv)/jT G Vh is uniquely determined by 


/ (divdiv);iT da: = u) V v GVh- 
Jq 

To present the commutative diagram, we need some interpolation operators. Let 
Qh be the orthogonal projection operator from L‘^{Vl) onto Vh which can be ex¬ 
tended to -G Vh as Vh C i7Q(H). 

For any element K G Tht define Ik ■ H'^{K) -G Pr{K) in the following way (cf. 
[5, 24, 22, 40]) : given w G any vertex a of K, and any edge e of K, 


{w 


Ikw{o) = w{o), 
lKw)'vds = 0 y V G Pr-2(e), 



Ikw)v da; 


0 Vug Pr-3{K). 
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The associated global interpolation operator Ih is given by 

iIh)\K := Ik for all K €%■ 


Let Ik = Ik Ik, Ih = Ih 'x h- 

Lemma 2.5. (divdiv)^ is a conforming discretization of B in the sense that 
ker((divdiv)?t) C kerS. 

Proof By the definition of It, we have (cf. [5, p. 1058]) 

(2.8) b{Th, v) = biru, hv), 'iru GVh,v G V. 

For any x G ker((divdiv)/j), we get from (2.8) that for any v G V, 


{Bt,v) = b{T,v) = b{T,Ifiv) = / (divdiv)/jX I/iri dx = 0. 

Jn 

Thus T G keiB. □ 

Then define Hk ■ H^{K,S) —)■ Pr-i{K,S) in the following way (cf. [5, 24, 22, 
16]): given r € H^{K,S), for any element K GTh and any edge e of K, 


Mn ((t - I1kt)\k) h-As 


0 VAiGPr-i(e), 


(t — TIkt) : ^dx = 0 


V <; G Pr-2(II,S). 


Jk 

The associated global interpolation operator lih : V —> V/i is given by 

{Tlh)\K:=IlK foralliFGT^. 

From the definition of II/i, it holds that 

(2.9) 5(t — II/iT, x) = 0 WTGV,vGVh- 

Namely QhB = (divdiv)/^!!/^. 

Lemma 2.6. The following sequence for the HHJ method 

_ C V®x (divdiv),! 

(2.10) Pi(L!;M")-- >Vh -^0 


is an exact sequence. 

Proof. As (2.4), (2.10) is a complex by direct computation. Then we prove 
ker((divdiv);i) = V® x Sh- Take any r G ker((divdiv)?i). Since ker((divdiv);i) C 
keri? and thus using (2.5) and the exact sequence (2.6) in the continuous level, we 
find a vector function cf) G II^(fl;R^) satisfying t = V® x 0. By direct computation, 
it hold for each K GTh 

curl(div(^lK)) = 2div(Tlif) G Pr-2{K,R?). 

Hence div(</)jif) G Pr-i{K), which combined with x cf) = t G Pr-i{K,S) means 
V(</>lif) G Pr-i{K,S). Therefore cf>\K G Pr{K,M.^), i.e. cf G Sh. 

Using the similar argument as in Lemma 2.1, we have ker(V® x) = Pi(U; K^). To 
show that (2.10) is exact, we shall prove (divdiv)/i(V;i) = Vh by adapting a technique 
in [5, p. 1056]. 






For any p GVh, let Wh & Vh be the solution of 

/ Vu>ft,-Vudx = — / pvdx \l V aVh- 


Let (To = 


. Thanks to M„((To) = ri^cron = Wh and Wh € Vh, ctq G V. 


Wh 0 
0 Wh 

Let (Tj = n^iCTo G Vh- Using (2.9), integration by parts twice, and the definitions of 
(Tg and Wh, it holds for any v G Vh 


b{(Ti,v) =b{(To,v) = ^ / (ToiV^udx- ^ / M„((To)5nuds 


K^Th 


= ^ / w/tAuds - ^ 


K^Th 
WhdnV ds 


K(^Th 


KeTh 


IdK 


= — / Vw/j-Vuda:= / pvdx, 

JQ Jq 

from which we can see that p = (divdiv)/icr/. The proof is finished. □ 

Theorem 2.7. We have the following commutative diagram for the HHJ method 


Pi{^; 


Pi(L!; 


Ih 

I 

-^Sh - 


ff^(divdiv,U;S) 


V“x 


-Vh 


(divdiv)h 


Qh 

-^Vh — 


■0 


0 


Proof. The identity Q?idivdiv = (divdiv);!!!^, has been proved in (2.9). 

Next we show that for any cj) G iT^(U;IR^) n dom(I?i), V® x {Ih4>) = LE/jV® x 4>. 
For each ^ G Pr- 2 {K,^) and K G Th, d follows from integration by parts and the 
definitions of Tlh and Ih 

(2.11) f X {Ihcf) -UhiV^ X cj,)) : <;dx = f V® x (/;,</)-(/.): ^ dx = 0. 

JK Jk 

On each e G £h{K), by the definition of 11^, it holds for any p G Pr-i{e) 

^ M„(V® X (Ihcj,) - n^(V® X cf,))p ds = ^ M„(V® X {Ihcj, - cf,))p ds. 

Note the fact that M„(V® x {Ih(t> — 4^)) = dtjjlhcf — <A) • n). Hence we get from 
integration by parts and the definition oi Ih 

(2.12) ^ M„(V® X jlhcj,) - n/,(V® X 4)))p ^s= dtiilhcj) - cf) ■ n)p ds = 0. 

Since (V® x (IhCp) — n?i(V® x (j)))\K £ Pk-i{K,^), (2.11)-(2.12) together with the 
wellposedness of 11^ means V® x {Ih4>) — n/i(V® x cj>) = 0, i.e. V® x {Ih4>) = 
n^(V® X cj,). □ 

Remark 2.8. It is worth to mention that we use the natural Sobolev spaces with 
minimal regularity in the top sequence of the commutative diagram. The interpolation 
operators Ih and II/j, however, are defined for smoother functions and not bounded in 
the corresponding Sobolev norms. Namely we treat these interpolation operators as 
densely defined unbounded operators. It is possible to use the smoothing procedure 
[3] to dehne stable quasi-interpolation operators while preserving the commutative 
property. □ 
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3. Stable Decomposition and Strengthened Cauchy Schwarz Inequality. 

In this section, we will present a stable decomposition for the space used in the 
HHJ method. We assume that there exists a sequence of meshes Ti,T 2 , ■ ■ ■ ,Tj =Th- 
Hereafter subscript k is used to indicate spaces associated to triangulation 7fc. The 
triangulation 7i is a shape regular triangulation of H and Tk+i is obtained by dividing 
each triangle in Tk into four congruent small triangles. The mesh size of Tk will be 
denoted by hk- By the construction, the ratio 7 ^ = hk+i/hk = 1/2. 

Based on the exact sequence (2.10), define JCk '■= V® x 5^ for A; = 1, 2, • • • , J. 
Obviously we have the following macro-decomposition 

}Ch — T ^2 “l- • • ■ -l- ACj. 

Denote by Nk the number of vertices in 7/ for k = 1, 2, • • • , J. Define the f-th 
patch uJk,i in the fc-th level as the union of the elements sharing the common f-th 
vertex in Tk for z = 1,2, • • • , Nk- Let 

Sk,i ■= G : supp(0) C ujk,i}^ ■= {t G Vfc : supp(x) C UJk,i}, 

and lCk,i ■= V® x Sk,i- It can be verified that 

./ J Nk J J Nk 

^ ^ ^ Sk,i, ^ Vfc = ^ ^ Vk,^, 

k—1 k — 1 i—1 k—1 k—1 i—1 


J J Nk 

( 3 . 1 ) ]Ch = Y,iCk = Y,Y.^'^N 

k—1 k—li—1 

We shall prove the space decomposition (3.1) is stable in the energy norm introduced 
by £^. 

3.1. Equivalent norms. We first introduce the following quotient spaces 


5 := I^G iT^(D;K2) . f ^da; = 0, f ^•a;da; = oj, 
I Jn Jn } 

:=<</) G : / cj)dx = 0, / cj) ■xdx = o\. 

I Jq Jn j 

It is easy to see that 

= SePi(fi;R^), Sk = Sk ePi(P,R^). 


Notation © means the direct sum. Since linear polynomial belongs to the space Sk, 
the spaces Sk are nested. Let 

W:={ct>eH\n-R^): [ cj)-ipdx = o w ip ePf°\p,R^)}, Wk-=NnSk. 

Jn 

It is obvious that cp^ G 5 if (/> G W, and vice versa. _ 

The following lemma says that in the quotient space 5, the differential operator 
V® X introduces a norm equivalent to norm. Similar result has been proved in [17] 
on a slightly different quotient space. 

Lemma 3.1. It holds 

ll0lli<l|V*x0|lo VcPgS. 


(3.2) 
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Proof. By a direct computation, we have for any vector cf and 


(3.3) 


V'* X 0 : V® X -0 = : e(t/;-^). 


Since 0 € 5, we have S W. According to the Korn’s inequality (see (2.2) in [12] 
and Theorem 2.3 in [21]), it follows 


ll‘/>^lli<lk(‘/>^)llo. 


Then we obtain from (3.3) 


ll</>lli = ll‘/>-"lli<l|£(0^)llo = llV*x0|lo, 


which ends the proof. □ 

3.2. Strengthened Cauchy Schwarz Inequality. Thanks to the relation (3.3), 
the following strengthened Cauchy Schwarz (SCS) inequality can be proved using the 
technique for the scalar case; see Xu [47]. 

Lemma 3.2. Let 1 < k < I < J. We have 


X 0 : V* X ^d:r < x </.|io||^||o V 0 S 5^,^ S Si. 


Next we prove the SCS inequality for the space decomposition (3.1) of K.. For 
this, we use the lexicographical order of the double index, i.e., (Z,j) > {k,i) if I > k 
or I = k, j > i. 

Theorem 3.3 (SCS). For any Tk,i € ICk.i and G JCij, we have 



Proof. Let Tk,i = V® x and = V’’ x ifij with 4>k,i ^ and -ipij G Sij. 
Using Lemma 3.2 and the fact that hf^\\ipij\\Q ~ [[V® x 'i/’iyllo; we get 



= 1 i=l l>k j=i 

./ Wfc Ni 


■J Wfc Ni 




fc=l i—l l>k j—1 
J Nk Ni 


^EEEEV-'^ii^^ X ‘/'mIIoIiv^ X 


fc=l i—l l>kj—1 
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On the other hand, since the index set nk(i) := {j G {* + Ij • ’ ’ > Nk},ujk,i Ciuikj ^ 0} 
is finite in the fcth level, 


J Nk Nk p J Nk p 

'■ ^k,j dx / Tk^i '■ <;k,j dx 

fc=l i=l jGnfc(i) ^ 


k—1 i—1 j=i-\-l ' 


J Nk 


1/2 


./ Ni 


1/2 


< 




\k—l i—1 


^1^1 j = l 

The summation of the last two inequalities implies the desired result. □ 

3.3. Stable Decomposition. Let Qf. be the projection from M^) onto 

Sk- It is easy to see that Q].4> S 5^ if 0 S 5. Due to the nestedness of spaces Sk, we 
also have QkQi — Qk for I > k. The following first order error estimate of Qk is well 
known 

(3.4) \\{I-Qkm<hkmi, ^oral\^P&H\n■,R^). 

Lemma 3.4. Let Wi = {Q^ — Qi_i)Sh for i = 1,2, • • • , J. ITe have 
[ V" X </> : V" X da; < x (t>\\o\\V" x xp\\o 

Jn 

for any cf) £Wi and G Wj. 

Proof According to the estimate of Qj-i and (3.2), 

mo = ii(j- g,_i)iMio < h,mi < hjim x i’h v i/. e w,. 

The proof is finished from Lemma 3.2. □ 

Let Pk be the V®x-orthogonal projection onto Sk, that is for any cf) G S, 


(3.5) 


V® X (Pfc0) : V" X xdx = / V" X 0 : V® X xda: V x G 5^. 


To derive the error estimate of Pk, we introduce another operator Rk which is related 
to the pure traction problem in_J;he planar linear elasticity. Let Rk ■ W —Wk be 
defined as follows: for any (p GW, Rk4> is uniquely determined by 

[ £{Rk4>) ■■ £{x) dx= [ £{cf)) : £{x) dx V x G Wk- 
Jn Jn 

According to the standard finite element approximation theory (cf. [9, (5.9)]), we have 

(3.6) ||0-Pfe</>||i_„</i^|10||i y<PGW 

for some constant a G (0,1]. Here a is the parameter indicating the elliptic regularity 
of the pure traction problem in the planar linear elasticity defined in Ll (cf. [27]). 
a = 1 if D is convex and 0 < a < 1 if 12 is nonconvex. 

Lemma 3.5. It holds 

\\cf>-PkcP\m<ht\\cP\\i VcPgS. 


(3.7) 
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Proof. Due to (3.3), (3.5) is equivalent to 


[ £{{Pk(t>)^) ■■ £{x^)dx = [ e{4>-^) : £{x^)dx V x € 5fc, 

Jci Jci 

which is nothing but 

[ £{{Pk(t))-^) ■■ e{x)dx = [ £{4>-^) : £{x)dx V x e Wfe- 
Jfi Jq 

Noting the fact that G W and {Pk4>)'^ G Wk, we get {Pk4>)'^ = Rki4>^). There¬ 
fore it follows from (3.6) 

\\<f- PkCfWl-a =||<A'^ - {Pk<f)^\\l-a = - i2fc(0‘^)||l-a 


as required. □ 

Again using the technique for the scalar space [47], we have the following 
stable decomposition of functions in Sh- 

Lemma 3.6 (Stable macro-decomposition). For each (p G St, there exists (pj. G 
Sk, k = 1,2, ■ ■ ■ ,J such that 


(P = '^4>k^ and ^ ||V" X ||V^ X 






Proof Let = Qk - Qfe-i. 4^k = Qk4> and xp^ = {P, - P^-l)p for i,k = 
1, 2, • • • ,J. Using Cauchy-Swarchz inequality, it holds 

k—l k—1 

,7 .7 

= E E / ^ X (Q,^,-)dx 

k—li,j — k ^ 

J iAi 

= E E / ^ : V* X (Q,^^.) dx 

i,j=i fc=i 
J iAj 

< E E 11^' (QfeTAJIIollV^ X iQk^P,)\\o, 

i,j — l k—l 


where i/\j = min{7, j}. According to the inverse inequality, the error estimate of Q^., 
and (3.7), we have 

l|v^ X {Qk^PMo < \Qk^r\i < hmQMi-c. < h-r\mW-o. < 

Combining the last two inequalities, we get from the strengthened Cauchy-Swarchz 
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inequality and (3.2) 


.7 ,7 iAj 

E11^^ X Ml < E E 

k—l 2J—1 k—1 

J J 

jj'=l i,i=l 

<E 11^*11? ^E11^'xV’.iig = iiv*xc/>iig. 

2—1 2—1 

On the other side, it follows from Lemma 3.4 and the strengthened Cauchy- 
Swarchz inequality 


IIV* X = E / V* X {Q,4>) : V* X (g,0)dx 
J 

<E7'*-''l|V^x(4c/,)||o||V^x(g,c/-)||o 

J .7 

<EI|V^x(Q,</.)||o = EI|V^x0J|o. 

z=l i=l 


The proof is completed. □ 

Lemma 3.7 (Stable micro-decomposition). Let (fif. = {Qf. — Qf._i)cj) with (p G Sh- 
Then based on the decomposition (3.1), there exists (pf, j S Sk,i, i = 1,2, ■ ■ ■ , Nk such 
that 


Nk Nk 

= Eii^'x0,,j|2<iiv*x0,i|2. 

i-1 i=l 

Proof. Let cpi. — <Pk,j be a decomposition such that supptpj^ j G ^k,j- Such 
decomposition can be obtained by partitioning the nodal basis decomposition of 4>k- 
For example, for a basis function associated to an edge, it can be split as half and 
half to the patch of each endpoint of this edge. 

According to the inverse inequality and the stability of the basis decomposition 
in L^-norm, we have 


Nk Nk 

E X cP^X ^ hf E ll‘^fc,*llo ^ h-X\Ml- 

2—1 2=1 

Since cpf. = {I — Qf._i)(pf., it holds from the estimate of Qk_i and (3.2) 
ll^fcllo = 11(1 - Qk-Mo < hkWcPdi < hkW^^ X cP,\\o. 

Therefore we can finish the proof from the last two inequalities. □ 

Hence the following multilevel stable decomposition of K.h can be derived by 
combination of Lemmas 3.6-3.7. 
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Theorem 3.8 (Stable decomposition). For each <t G K-h, there exists crk,i G JCk,i, 
fc = 1, 2, • • • , J, i = 1, 2, • • • , N^. such that a 

J Nk J Nk 

^ = and - ll'^llo- 

k—1 i—1 k—1 i—1 

Proof. Since cr G /C/i, we can find a unique element cf G Sh such that <j = V® x (p. 
Let (pi^ = {Qi^ — Qj^_i)(p. We get from Lemma 3.6 

,/ J 

(3.8) = and ^ ||V^ x ||V* x 0||g. 

k=l k=l 

Then we apply Lemma 3.7 to obtain a decomposition of (pf. such that 
Nk W 

(3.9) 5]||V«x</.fc,J|g<||V*x0,||2 

i—1 i—1 

with (pi^ j e Sk^i for z = 1, 2, • • • , iV^ and fc = 1,2, • • • , J. Now let cr^, ^ = V® x ^ S 
ICk.i- It is apparent that 


J Nk 

O’ = ^ ^ (Tk^i. 

k—1i—1 

Moreover, it follows from (3.8)-(3.9) 

./ Vfc ,7 ATfe J 

E E = E E11^' < E ^ X ‘^iio = ikiio- 

k—1 i—1 k—1 i—1 k—1 

The proof is ended. □ 

4. Multigrid Methods for the HHJ mixed methods. In this section we 
shall develop a multigrid method using an overlapping Schwarz smoother for the 
HHJ mixed method and prove its uniform convergence. We first solve a Poisson 
equation with Dirichlet boundary condition to transfer the source. Then we apply 
the multilevel method advised in [19] and the space decomposition (3.1) to obtain a 
V-cycle multigrid method with an overlapping Schwarz smoother for the HHJ mixed 
method. We analyze the V-cycle multigrid method by using the stable decomposition 
and the strengthened Cauchy Schwarz inequality. 

4.1. Reformulation. We change the source to the first equation in the saddle 
point system (2.2)-(2.3). One possibility is as follows: let Wh G Ph be the solution of 


Vwh ■ Vvh dx= fvh da; ^ Vh G Th¬ 


in 


In 


This is the standard Poisson equation which can be solved efficiently by multigrid 

Wh 0 


methods. Let erg = 


According to the proof of Lemma 2.6, we have 


0 Wh 

Mn{cro) = Wh, (Tg € V, UhfTo G Vh and (divdiv);,n/,crg = -Qhf, i.e., 


b{Ilhcro,Vh) = - fvhdx \/vhGVh- 


n 
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Now set (Th — (Th + Il/icro, then the HHJ mixed method (2.2)-(2.3) is equivalent to: 
Find (S-fiyUh) GVh ^-Vh such that 

(4.1) a{a-h,r) + b{T,Uh) = -a(Uhcro,T) 

(4.2) b{a-h,v) =0 VvGVh- 

After obtaining (Th, due to Theorems 5.1-5.2 in [32], we can acquire deflection by 
solving the following Poisson problem using standard multigrid methods: Find Uh G 
Vh such that 

/ Vuh ■ Vvhdx = a{(Th,'nhTo) Mvh^Vh- 

Jn 

"“■""=( 0 .")■ 

Our multigrid method is actually developed for solving (4.1)-(4.2). 

4.2. A V-cycle Multigrid Method. We shall use the multilevel methods for 
constrained minimization problems developed in [19] and adapt to the HHJ mixed 
method under consideration. For simplicity, we consider the lowest order HHJ method 
for which Vh consists of piecewise constant symmetric matrix function and normal- 
normal component is continuous, Sh is the standard linear finite element space for 
vector functions, and Vh is the linear finite element space with zero boundary condi¬ 
tion for scalar functions. For high order HHJ methods, we can combine the multigrid 
cycles for the lowest order and an overlapping Schwarz smoother in the finest level to 
design efficient multigrid solvers. 

Let Mk : Vk —>■ Vk be the mass operator associated with the bilinear form a{-, •): 
for any t €Vk, MkT G Vk is uniquely determined by 


/ MfcT : ^dx = a(x, <;) V ^ S Vfc. 

Jn 

The mixed variational problem in the fc-th level is: Find {o'k, Uk) G Vk x Vk such that 


a{(Tk,T) + b{T,Uk) = rk-rdx VrSVfc, 

(4.3) Jn 

b{a-k,v) =0 VvGVk, 

with the residual rk G L‘^(fi;S). 

As we mentioned before, the smoother in each level is an overlapping Schwarz 
method. To simplify the notation, we skip the level index k and describe the local 
problem in each subspace Vi (of a given level k) below. Define Mi : Vi ^ Vi as 
for (Ti G Vi, Mrxi G Vi such that (Micri,Ti) = {M(Ti,Ti) for all Ti G Vi- Let 
Vi = V C\ divdiv(Vi). Define Bi : Vi -G Vi as for cTi G Vi, divdivcTi G Vi such that 
{Bi(Ti,Vi) = (divdiv(Tj,Ui) for all Vi G Vi. 

(«) (“: ?)(::) 

Let uJi be the support of Vi- For the lowest order HHJ method, this is the patch 
of the i-th vertex of the triangulation in the given level. The space Vi is spanned 
by basis functions associated to all edges connecting to the i-th vertex. The matrix 
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representation of Mi can be extracted from the global one using the edge index in Wi. 
The right-hand side is the corresponding components of / minus the contribution from 
the current approximation. Note that MtTi-i only need to be computed locally by 
including the boundary edge index of duji. The correct space Vi is somehow difficulty 
to identify. We shall work on the space ICi instead. An algebraic way to find K.i is 
as follows. We extract a sub-matrix of Bi consisting of all nonzero entries associated 
to the edge index in Vi and compute ker(i?i) numerically. An alternative way is 
computing e^{<pi) where (pi is the vector hat function associated to vertex i. 

Remark 4.1. Since ker^Bh) = x Sh due to the exact sequence (2.10), the 
mixed method (4.1)-(4.2) can be rewritten as: Find cpf^ G such that 

(4.5) a(V'* X 0^, V* X ^) = -a{nh(To, V* x -0) 

with ffh = V® X By the theory in [34], this symmetric and positive semidefinite 
problem can be solved by multigrid method efficiently. Solving the mixed method 
(4.1)-(4.2) is essentially equivalent to the multigrid method developed for (4.5). □ 

We then discuss the prolongation and restriction operators. Since both finite 
element spaces Vk and Vk are nested, the prolongations Ik-i ■ Vk-i —t Vk and I^_i : 
Vk-i —t Vk are chosen as the natural inclusions. Set the restriction := 
and := . With the restriction and prolongation matrix, the matrices Mk 

and Bk in each level can be obtained by the standard triple product. 

With previous preparation, a V-cycle multigrid method for problem (4.1)-(4.2) is 
summarized in Algorithm 1 with rj = —MjTIj(Tq. 

The A-norm introduced by a(',') on Vh is equivalent to the L^-norm. With 
the stable decomposition and the strengthened Cauchy-Schwarz inequality proved in 
Section 3, applying the framework developed in [19], we concluded that multigrid 
method Algorithm 1 is a contraction with contraction number bounded away from 
one uniformly with respect to mesh size as follows. 

Theorem 4.2. Let [ahTUh) be the solution of the mixed method (4.1)-(4.2). 
Given an initial guess ctq € Vh, l^i cr be the kth iteration in Algorithm 1. Then there 
exists a constant 6 G (0,1) independent of the mesh level such that 

with \\t\W := a{T, t). 

4.3. Numerical Results. To confirm the theoretical results established in the 
previous sections, numerical experiments are carried out. The simulation is imple¬ 
mented using the MATLAB software package iFEM [18]. Starting from an initial 
grid, several uniform refinement are applied to obtain a fine mesh. The level listed 
in the first column indicates how many refinements applied and the size of the saddle 
point system is listed in the second column. The stopping criterion is the relative 
residual is less than 10“®. The iteration steps are reported in Table 4.1. 

We test two examples. One is a square SI = (0,1) x (0,1) and another is an 
L-shaped domain SI = (—1,1) x (—1,1)\[0,1) x (—1,0]. For the square domain, the 
Poisson ratio is ^ = 0.3 and the exact solution of (2.1) is chosen as 

u{x,y) = {x^ - xf{y'^ - yf. 

And for L-shaped domain, we simply set / = 1 and the Poisson ratio v = Q. The 
later example is to test the multigrid method for problems without full regularity 
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Algorithm: MG(A:, cr^, r^) 
if A: = 1 then 

solve problem (4.3) exactly; 

end 

if /c > 1 then 
Presmoothing 
for j = 1 : TOi do 
<7’fc,0 <^k] 

for i = \ ■. Nk do 

Update by solving local problem (4.4); 

end 

cTfc ■<— erk.Nk j 

end 

Coarse grid correction 

rk-i ^ - MfcCTfe); 

ef_i ^MG(fc-l,0,rfc_i); 

o’fc ^ cTfc + 

Postsmoothing 
for j — 1 : m 2 do 
cr k,0 CTk', 

for i = Nk : — 1 : 1 do 

Update a-k^i by solving local problem (4.4); 

end 

CTfc CTk^Nk j 
end 
end 

Algorithm 1: A V-cycle multigrid method for problem (4.1)-(4.2). 

assumption. From Table 4.1 we can see that the iteration steps of V-cycle multigrid 
method almostly remain invariable when the mesh size becomes smaller and smaller, 
as Theorem 4.2 indicates. Moreover through the comparison of different number of 
smoothings, we conclude that one smoothing is enough. Two smoothing steps will 
save only few iteration steps but with more computational cost since the cost of one 
V-cycle with 2 smoothing steps is almost doubled that with 1 smoothing step. This 
indeed shows the advantage of removing the assumption of requiring large enough 
smoothing steps. These numerical results are all in coincide with the theoretical 
result Theorem 4.2. 

5. Conclusion. In this paper, we have advanced and analyzed a V-cycle multi¬ 
grid method with an overlapping Schwarz smoother for HHJ mixed method. The 
novelties of our V-cycle multigrid method are: 

(1) Full regularity assumption is not necessary for our multigrid method, i.e. our 
approach works for both convex and non-convex domains. 

(2) One smoothing is enough to guarantee the uniform convergence of our V-cycle 
multigrid algorithm, whereas large enough smoothing steps are usually required 
in the former multigrid methods for the fourth order partial differential equation. 
To obtain the uniform convergence of our V-cycle multigrid algorithm, we estab¬ 
lish the exact sequence for the HHJ mixed method in both the continuous and discrete 
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Table 4.1 

Iteration steps of V-cycle multigrid for the saddle point system with {mi,m 2 ): mi pre¬ 
smoothing and m 2 post-smoothing steps. Stopping criterion is the relative residual is less than 
10“®. The left table is on the unit square example with i/ = 0.3 and the right one is the L-shaped 
domain example with u = 0. 


level 

size 

(1,1) 

(2,2) 

level 

size 

(1,1) 

(2,2) 

3 

1,089 

18 

14 

3 

833 

13 

11 

4 

4,225 

21 

15 

4 

3,201 

17 

14 

5 

16,641 

22 

16 

5 

12,545 

19 

16 

6 

66,049 

23 

16 

6 

49,665 

20 

17 


levels, and prove the stable decomposition and strengthened Cauchy Schwarz inequal¬ 
ity. Then using the framework developed in [19] we obtain the uniform convergence. 
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